2022-09-30 2D Advection-diffusion and waves
Contents
2022-09-30 2D Advection-diffusion and waves#
Last time#
Implicit Runge-Kutta methods
Exploring/discussing tradeoffs
SciML benchmarks suite (
DifferentialEquations.jl)
Today#
FD methods in 2D
Cost profile
The need for fast algebraic solvers
Wave equation and Hamiltonians
Symplectic integrators
using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays
function my_spy(A)
cmax = norm(vec(A), Inf)
s = max(1, ceil(120 / size(A, 1)))
spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end
function plot_stability(Rz, title; xlims=(-2, 2), ylims=(-2, 2))
x = LinRange(xlims[1], xlims[2], 100)
y = LinRange(ylims[1], ylims[2], 100)
heatmap(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2), aspect_ratio=:equal, title=title)
end
struct RKTable
A::Matrix
b::Vector
c::Vector
function RKTable(A, b)
s = length(b)
A = reshape(A, s, s)
c = vec(sum(A, dims=2))
new(A, b, c)
end
end
function rk_stability(z, rk)
s = length(rk.b)
1 + z * rk.b' * ((I - z*rk.A) \ ones(s))
end
rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0], [1, 2, 2, 1] / 6)
function ode_rk_explicit(f, u0; tfinal=1., h=0.1, table=rk4)
u = copy(u0)
t = 0.
n, s = length(u), length(table.c)
fY = zeros(n, s)
thist = [t]
uhist = [u0]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
for i in 1:s
ti = t + h * table.c[i]
Yi = u + h * sum(fY[:,1:i-1] * table.A[i,1:i-1], dims=2)
fY[:,i] = f(ti, Yi)
end
u += h * fY * table.b
t = tnext
push!(thist, t)
push!(uhist, u)
end
thist, hcat(uhist...)
end
ode_rk_explicit (generic function with 1 method)
Extending advection-diffusion to 2D#
1 dimension#
Cell Peclet number \(\mathrm{Pe}_h = \frac{\lvert w \rvert h}{\kappa}\)
\(\mathrm{Pe}_h \lesssim 1\) avoids oscillations
\(\mathrm{Pe}_h \gtrsim 1\) is non-stiff for time-dependent model
Centered versus upwind for advection
Need uniformly bounded \(\kappa \ge \epsilon > 0\)
“Strong form” not defined at discontinuities in \(\kappa\)
Works okay using divergence form and fluxes at staggered points
2 dimensions#
\(\Omega\) is some well-connected open set (we will assume simply connected) and the Dirichlet boundary \(\Gamma_D \subset \partial \Omega\) is nonempty.
Finite difference methods don’t have an elegant/flexible way to specify boundaries
We’ll choose \(\Omega = (-1, 1) \times (-1, 1)\)
On finite difference grids#
Non-uniform grids can mesh “special” domains
Rare in 3D; overset grids, immersed boundary methods
Concept of staggering is complicated/ambiguous

Time-independent advection-diffusion#
Advection#
If we choose divergence-free \(\mathbf w\), we can use the stencil
Diffusion#
When would you trust this decomposition?
If we have constant \(\kappa\), we can write
\[\begin{split} -\kappa \nabla\cdot \nabla u \approx \frac{\kappa}{h^2} \begin{bmatrix} & -1 & \\ -1 & 4 & -1 \\ & -1 & \end{bmatrix} \!:\! \begin{bmatrix} u_{i-1, j+1} & u_{i, j+1} & u_{i+1,j+1} \\ u_{i-1, j} & u_{i, j} & u_{i+1,j} \\ u_{i-1, j-1} & u_{i, j-1} & u_{i+1,j-1} \\ \end{bmatrix} \end{split}\]
Advection-diffusion in code#
function advdiff_matrix(n; kappa=1, wind=[1, 1]/sqrt(2))
h = 2 / n
rows = Vector{Int64}()
cols = Vector{Int64}()
vals = Vector{Float64}()
idx(i, j) = (i-1)*n + j
stencil_advect = [-wind[1], -wind[2], 0, wind[1], wind[2]] / h
stencil_diffuse = [-1, -1, 4, -1, -1] * kappa / h^2
for i in 1:n
for j in 1:n
if i in [1, n] || j in [1, n]
push!(rows, idx(i, j))
push!(cols, idx(i, j))
push!(vals, 1.)
else
append!(rows, repeat([idx(i,j)], 5))
append!(cols, [idx(i-1,j), idx(i,j-1), idx(i,j), idx(i+1,j), idx(i,j+1)])
append!(vals, stencil_advect + stencil_diffuse)
end
end
end
sparse(rows, cols, vals)
end
advdiff_matrix (generic function with 1 method)
Spy the matrix#
A = advdiff_matrix(10, wind=[1, 3], kappa=.1)
my_spy(A)
A = advdiff_matrix(20, wind=[0, 0], kappa=.1)
ev = eigvals(Matrix(A))
scatter(real(ev), imag(ev))
Plot a solution#
n = 100
x = LinRange(-1, 1, n)
y = x
f = cos.(pi*x/2) * cos.(pi*y/2)'
heatmap(x, y, f, aspect_ratio=:equal)
A = advdiff_matrix(n, wind=[3,1], kappa=.01)
u = A \ vec(f)
heatmap(x, y, reshape(u, n, n), aspect_ratio=:equal)
What happens when advection dominates?
As you refine the grid?
Cost breadown and optimization#
using ProfileSVG
function assemble_and_solve(n)
A = advdiff_matrix(n)
x = LinRange(-1, 1, n)
f = cos.(pi*x/2) * cos.(pi*x/2)'
u = A \ vec(f)
end
@profview assemble_and_solve(200)
┌ Warning: The depth of this graph is 98, exceeding the `maxdepth` (=50).
│ The deeper frames will be truncated.
└ @ ProfileSVG /home/jed/.julia/packages/ProfileSVG/ecSyU/src/ProfileSVG.jl:275
┌ Warning: The depth of this graph is 98, exceeding the `maxdepth` (=50).
│ The deeper frames will be truncated.
└ @ ProfileSVG /home/jed/.julia/packages/ProfileSVG/ecSyU/src/ProfileSVG.jl:275
What’s left?#
Symmetric Dirichlet boundary conditions
Symmetric Neumann boundary conditions
Verification with method of manufactured solutions
Non-uniform grids
Upwinding for advection-dominated problems
Variable coefficients
Time-dependent problems
Fast algebraic solvers
Gas equations of state#
There are many ways to describe a gas
Name |
variable |
units |
|---|---|---|
pressure |
\(p\) |
force/area |
density |
\(\rho\) |
mass/volume |
temperature |
\(T\) |
Kelvin |
(specific) internal energy |
\(e\) |
energy/mass |
entropy |
\(s\) |
energy/Kelvin |
Equation of state#
Ideal gas#
pressure(rho, T) = rho*T
contour(LinRange(0, 2, 30), LinRange(0, 2, 30), pressure, xlabel="\$\\rho\$", ylabel="\$T\$")
Conservation equations#
Mass#
Let \(\mathbf u\) be the fluid velocity. The mass flux (mass/time) moving through an area \(A\) is
If mass is conserved in a volume \(V\) with surface \(A\), then the total mass inside the volume must evolve as
where we have applied the divergence theorem. Dropping the integrals over arbitrary volumes, we have the evolution equation for conservation of mass.
Momentum#
The momentum \(\rho \mathbf u\) has a flux that includes
convection \(\rho \mathbf u \otimes \mathbf u\)
this is saying that each component of momentum is carried along in the vector velocity field
pressure \(p I\)
viscous \(-\boldsymbol\tau\)
A similar integral principle leads to the momentum equation
Simplifications#
Ignore viscous stress tensor \(\boldsymbol \tau\)
Ignore energy equation (not yet written) and assume constant temperature
\(p = a^2 \rho\) where \(a\) is the acoustic wave speed
Linearization#
Split each state variable into a mean state and a small fluctuation
\(\rho = \bar\rho + \tilde\rho\)
\(u = \bar u + \tilde u\)
Let \(\widetilde{\rho u} = (\bar\rho + \tilde\rho) (\bar u + \tilde u) - \bar\rho\bar u \approx \tilde \rho \bar u + \bar\rho \tilde u\), where we have dropped the second order term \(\tilde \rho\tilde u\) because both are assumed small.
We consider background state \(\bar u = 0\) and constant \(\bar\rho(x,y,t) = \bar\rho\). Then
Two forms of acoustic wave equation#
Divide the momentum equation through by background density and dropping the tildes yields the standard form.
Examine second equation
Let’s differentiate the first equation,
Note: we had to assume these derivatives exist!
We can reduce this to a first order system as
Question#
How is the problem size different?
What might we be concerned about in choosing the second formulation?
Laplacian in periodic domain#
function laplacian_matrix(n)
h = 2 / n
rows = Vector{Int64}()
cols = Vector{Int64}()
vals = Vector{Float64}()
wrap(i) = (i + n - 1) % n + 1
idx(i, j) = (wrap(i)-1)*n + wrap(j)
stencil_diffuse = [-1, -1, 4, -1, -1] / h^2
for i in 1:n
for j in 1:n
append!(rows, repeat([idx(i,j)], 5))
append!(cols, [idx(i-1,j), idx(i,j-1), idx(i,j), idx(i+1,j), idx(i,j+1)])
append!(vals, stencil_diffuse)
end
end
sparse(rows, cols, vals)
end
laplacian_matrix (generic function with 1 method)
L = laplacian_matrix(10)
ev = eigvals(Matrix(L))
scatter(real(ev), imag(ev))
Wave operator#
function wave_matrix(n; a=1)
Z = spzeros(n^2, n^2)
L = laplacian_matrix(n)
[Z I; -a^2*L Z]
end
wave_matrix(2)
8×8 SparseMatrixCSC{Float64, Int64} with 16 stored entries:
⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0
-4.0 2.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅
2.0 -4.0 ⋅ 2.0 ⋅ ⋅ ⋅ ⋅
2.0 ⋅ -4.0 2.0 ⋅ ⋅ ⋅ ⋅
⋅ 2.0 2.0 -4.0 ⋅ ⋅ ⋅ ⋅
A = wave_matrix(8; a=2) * .1
ev = eigvals(Matrix(A))
plot_stability(z -> rk_stability(z, rk4), "RK4", xlims=(-4, 4), ylims=(-4, 4))
scatter!(real(ev), imag(ev), color=:black)
Question: would forward Euler work?#
Example 2D wave solver with RK4#
n = 20
A = wave_matrix(n)
x = LinRange(-1, 1, n+1)[1:end-1]
y = x
rho0 = vec(exp.(-9*((x .+ 1e-4).^2 .+ y'.^2)))
sol0 = vcat(rho0, zero(rho0))
thist, solhist = ode_rk_explicit((t, sol) -> A * sol, sol0, h=.02)
size(solhist)
(800, 51)
@gif for tstep in 1:length(thist)
rho = solhist[1:n^2, tstep]
contour(x, y, reshape(rho, n, n), title="\$ t = $(thist[tstep])\$")
end
┌ Info: Saved animation to
│ fn = /home/jed/cu/numpde/slides/tmp.gif
└ @ Plots /home/jed/.julia/packages/Plots/lW9ll/src/animation.jl:137
Accuracy, conservation of mass with RK4#
thist, solhist = ode_rk_explicit((t, sol) -> A * sol, sol0, h=.05,
tfinal=1)
tfinal = thist[end]
M = exp(Matrix(A*tfinal))
sol_exact = M * sol0
sol_final = solhist[:, end]
norm(sol_final - sol_exact)
0.020151117484358525
mass = vec(sum(solhist[1:n^2, :], dims=1))
plot(thist[2:end], mass[2:end] .- mass[1])
Conservation of energy with RK4#
Hamiltonian structure#
We can express the total energy for our system as a sum of kinetic and potential energy.
where we identify \(\rho\) as a generalized position and \(\dot\rho\) as generalized momentum. Hamilton’s equations state that the equations of motion are
where we have used the weak form to associate \(\int \nabla v \cdot \nabla u = v^T L u\).
function energy(sol, n)
L = laplacian_matrix(n)
rho = sol[1:end÷2]
rhodot = sol[end÷2+1:end]
kinetic = .5 * norm(rhodot)^2
potential = .5 * rho' * L * rho
kinetic + potential
end
ehist = [energy(solhist[:,i], n) for i in 1:length(thist)]
plot(thist, ehist)
Velocity Verlet integrator#
function wave_verlet(n, u0; tfinal=1., h=0.1)
L = laplacian_matrix(n)
u = copy(u0)
t = 0.
thist = [t]
uhist = [u0]
irho = 1:n^2
irhodot = n^2+1:2*n^2
accel = -L * u[irho]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
u[irho] += h * u[irhodot] + h^2/2 * accel
accel_next = -L * u[irho]
u[irhodot] += h/2 * (accel + accel_next)
accel = accel_next
t = tnext
push!(thist, t)
push!(uhist, copy(u))
end
thist, hcat(uhist...)
end
wave_verlet (generic function with 1 method)
thist, solhist = wave_verlet(n, sol0, h=.04)
@gif for tstep in 1:length(thist)
rho = solhist[1:n^2, tstep]
contour(x, y, reshape(rho, n, n), title="\$ t = $(thist[tstep])\$")
end
┌ Info: Saved animation to
│ fn = /home/jed/cu/numpde/slides/tmp.gif
└ @ Plots /home/jed/.julia/packages/Plots/lW9ll/src/animation.jl:137
Accuracy and conservation for Verlet#
thist, solhist = wave_verlet(n, sol0, h=.05, tfinal=50)
tfinal = thist[end]
M = exp(Matrix(A*tfinal))
sol_exact = M * sol0
sol_final = solhist[:, end]
@show norm(sol_final - sol_exact)
mass = vec(sum(solhist[1:n^2, :], dims=1))
plot(thist[2:end], mass[2:end] .- mass[1])
norm(sol_final - sol_exact) = 6.86250099252763
ehist = [energy(solhist[:,i], n) for i in 1:length(thist)]
plot(thist, ehist)
Notes on time integrators#
We need stability on the imaginary axis for our discretization (and the physical system)
If the model is dissipative (e.g., we didn’t make the zero-viscosity assumption), then we need stability in the left half plane.
The split form \(\rho, \rho\mathbf u\) form is usually used with (nonlinear) upwinding, and thus will have dissipation.
Runge-Kutta methods#
Easy to use, stability region designed for spatial discretization
Energy drift generally present
Verlet/leapfrog/Newmark and symplectic integrators#
These preserve the “geometry of the Hamiltonian”
energy is not exactly conserved, but it doesn’t drift over time
such methods are called “symplectic integrators”
May not have stability away from the imaginary axis (for dissipation)
Most require a generalized position/momentum split, “canonical variables”